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Abstract 

The behavior of nano-voids composed of vacancies (Vs) at grain boundaries (GBs) is fundamental to 
the design of the radiation tolerance of poly-crystalline metals (PCs) via GB engineering. In this study, based 
on differential evolution, a framework for determining the stable structure of GB nano-voids is developed. 
Combining the framework with multiscale simulations, we elucidate the vacancy-accumulation and GB void 
formation mechanism under irradiation. A GB-structure dependent picture is revealed. At special 
coincidence-site-lattice (CSL) GBs of 25(310) and 25(210) with a medium V-GB binding energy, the V 
could be reemitted from the GB and also has driving force to be clustered at the GB, developing particularly 
stable V-clusters from a linear configuration to a platelet and finally to three-dimensional void that has large 
strain fields in iron with small bulk modulus and a bulk-void alike structure in the GB with large bulk 
modulus. A group of vacancies reconstruct their positions during the growth. The ripening is also mediated 
by the mobility of small V-clusters in addition to free Vs. General high-angle and low-angle GBs trap Vs 
efficiently, where V-clusters only align one-dimensionally or hardly nucleate. Based on the bonding among 
the vacancies and their neighboring atoms of a nano-void, we propose a high-accuracy predictive linear 
energetic model applied to the nano-void both at the iron/molybdenum/tungsten GBs and in the grain 


interior. The model captures the anisotropic feature of a nano-void and reproduces the oscillated vacancy 


energy level near a nano-void, showing distinct advantages over conventional continuum model and Wulff 
construction based energy model. Finally, the collective behavior of multiple GBs plays a role in the GB 
void formation. The present work offers fundamental mechanistic insights to GB nano-void formation and 


growth and sets a key step towards GB-void prevention in PCs by reducing the fraction of special CSL-GBs. 
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1. Introduction 

Nano-voids composed of vacancies (Vs) and their clusters (Vas), one of natural native defects prevalent 
in materials under radiation [1—4], deformation [5], annealing [6], and high hydrogen pressure [7], play an 
important role in determining materials properties and performance, e.g. swelling [1—4,8,9], creep [10,11], 
impurity retention [12—15], and intergranular corrosion [16]. Notably, the interface e.g. grain boundaries 
(GBs) serving as sinks for defects appears to act as nucleation sites for vacancies to coalesce. The voids (or 
cavities containing a limited amount of gaseous impurities) have been observed to form preferentially along 
GBs and at triple junctions in poly-crystalline metals (PCs) [2-4,10—14,16,17]. Since the GB voids observed 
as early as 1954 [17], the thermodynamic models for vacancy accumulation and void growth at GBs have 
been proposed [18,19]. However, the fundamental properties for the GB nano-void (e.g. its structure, 
energetics and kinetics) are still lacked except the properties of small V,s in bulk metals [20-25], which 
necessitate thorough studies on the behavior of GB nano-voids. 

Besides, the GB structure probably affects the nano-void formation mechanism in PCs [2-4,26,27] by 
providing different vacancy sink strengths and in-boundary properties of a nano-void. Experiments of the 
neutron and electron irradiations to PC Fe-15Cr-15Ni steel and PNC316 stainless steel at 749-775 K suggest 
that the void-denuded zone (VDZ) was formed only near a random GB, while voids were formed both in the 
grain interior and at coincidence site lattice (CSL) GBs [3]. In electron irradiation experiments at 573 K on 
Cu thin films [2], very few voids were observed around high-angle random GBs and low-angle GBs, while 
voids were preferentially observed to nucleate and grow in the regions with a high density of the coherent 
twin boundaries. The difference has been speculated to arise from the low sink strength of the low energy of 
the CSL-GBs [26]. On the other hand, however, the radiation response of NCs with high-density of GBs 
exhibits a scattered feature with respect to the GB character. For instance, VDZs were recently found to be 
independent of grain size, grain orientation, and GB misorientation angle in NC Fe [27]. Although the void 


size observed in these experiments is as large as tens of nanometers, the different large void formation 


conditions are possibly related to the early period of the nano-void formation at the GBs with distinct 
structures. Therefore, the GB character complicates the GB nano-void behavior. 

To obtain a full picture about a GB nano-void behavior, the first step is to find its stable structure, which 
acts as the basis of both the energetic and kinetic properties of the cluster. For this purpose, it in principle 
requires exploring the potential energy surface (PES) on the complex 3N-dimensional (N is the number of 
atoms in the system) configuration space and locating the global minimum of the PES. The challenging 
difficult exists as the cluster size is over a certain limit (e.g. 10 inaccessible to the enumerative algorithm) 
since the possible number of configurations increases exponentially with the cluster size. There already exist 
several methods for PES exploring, e.g. enumeration search methods [20—22,28,29] and the autonomous 
basin climbing method [30]. However, these methods may not be efficient enough to find the global 
minimum of the PES due to the exceptionally large configurational space of a Vy. Alternatively, in recent 
years, the biological intelligent global optimization algorithm e.g., genetic algorithm (GA), particle swarm 
optimization (PSO) and differential evolution (DE) have been proved to be successful in optimizing atomic 
clusters structure [31], crystal structure [32,33], interface/boundary structure [34] in the material science 
community and solving complex optimization problems in other fields [35]. The GA has been used to 
investigate the structure and stability of defects clusters in quite limited systems [36-38]. Compared to other 
group intelligence methods, the DE algorithm as one of heuristic group intelligence approaches to solve 
global optimization problems has advantages of simplicity, robustness and high performance [35,39]. 

To address the aforementioned challenge, in the present work, an approach in the framework of the DE 
was developed to find the stable configuration of a nano-void at Fe/Mo/W GBs and also in their grain 
interiors. The choice of the three model systems is based on the consideration of their identical crystal 
structure in addition to the nuclear industrial application of Fe/Mo/W as structural materials in nuclear 
reactors [26,40,41]. Besides, we also tried to relate the formation energy of a nano-void to its geometrical 


parameters (e.g. the bonding between V-V and the bonding between the V and the surrounding atoms). And 


then, a predictive discrete model for the GB nano-void energetics was built. In addition, the kinetics of a Va 
in the vicinity of the GB was investigated by using lattice kinetic Monte Carlo (LKMC) method [42]. 
Contrary to the object kinetic Monte Carlo (OKMC) with a pre-constructed rate catalog for the events [42], 
the rate catalog was rebuilt at every step of the LKMC. Furthermore, the energy barriers of all the transition 


paths of a nano-void were calculated in a parallel way by using nudge elastic band (NEB) method [43]. 


2. Simulation methods 


2.1. Atomic GB models 

The interatomic interaction between Fe atoms was described by the empirical potential proposed by 
Mendelev et al. as part of their embedded-atom method (EAM) (Potential 2 in Ref. [44]). The EAM 
potentials were also used to describe interatomic interaction in Mo [45] and W [46] (EAM4 in Ref. [46]). 
Since the vacancy behavior at the GB relates to the GB structure, multiple GBs were modeled in this work, 
including 44 symmetric tilt GBs of [A k 0]/[0 0 1] with tilt axis of [0 0 1] ( index A and k varies within the 
range of [2 13] and [1 11], respectively). The misorientation angle ( @) for these GBs varies from 8.8 to 84.6 
°; X-value varies from 5 to 265. For comparison, several GBs with another tilt axis of [1 1 0] were also 
modeled, including £3(1 1 1)/[1 1 0], £3(1 1 2)/[1 1 0], E9(1 1 4)/[1 1 0], 92 2 1)/[1 1 0], 2110 1 3)/f1 1 


0], and £11(3 3 2)/[1 1 0]. 


2.2. Details for searching the ground state configuration of a V, at the GB by using the differential evolution 
(DE) algorithm 

The life group intelligent algorithm DE [35,39] was employed to solve the global optimization issue 
related to finding the stable configuration of nano-void (Fig. 1(a)). The system was divided into two regions: 
the core region Regl where the V, locates and its environment Reg2 (Fig. 1(b)). The space near a Vn in 
Reg1 was divided into mesh. Each site of the mesh denotes a state of the V composing a V, in the un-relaxed 


space. Note that, the real structure of a Vn should be relaxed. Such notation of the V, maps the relaxed space 


to the un-relaxed space: 


unrelaxed (V) <> relaxed (V), R R +E, (1) 


relaxed ~~ *unrelaxed 


where R iretacca AA Roueg I the Vn configuration before and after relaxation, respectively; ¢ is 
relaxation-induced displacement of a Va. The individual V in a Vn before relaxation was considered to locate 
at a site near a lattice in the un-relaxed space (Eq. 1). For a Vn, the V composing the V, was considered to 
locate at a crystal lattice site. Particularly, for the small V, at the GB, the space mesh was restricted to the 
position sites of a single V in the V, (the site with minimum V formation energy). Such restriction could 
significantly improve the search efficiency. The above approximation was found to be rational by comparing 
the configurations of the small V, found via the DE method with that found by enumeration search method. 
Furthermore, multiple runs of the optimizer lead to the same stable configuration of a nano-void. For the 
large Vn, considering the possible expansion of the occupation range for the V in the V, towards the bulk 
region from the GB, the mesh was also restricted to the region near the GB rather than only within the GB 
plane. 

Based on above approximation, the problem of finding the global energy minimum on the PES of the Va 
system containing N/ sites in Regl was transformed into to a combinatorial problem of choosing n sites 
from N/ and minimizing the system energy: 


min E = E(x), x =[x1,x2,..., xn], xi | e[l, N1], (2) 


i=1,n 
where xi is an integer denoting the individual site label in Regl. The DE was employed to solve this 
optimization problem over discrete integer space, including four steps of initialization, mutation, crossover, 
and selection, as illustrated in Fig. 1(c). 
(a). Initialization 

For a given calculated radius r of Regl, the candidate sites set C={C1, C,..., Cui} in Regl. Then, the 


basic parameters of the DE were set, including the population size Np, scale factor Fscale, and crossover 


probability Per. For each individual x’ |-1.2,-.\p» It was assigned to n different integers randomly chosen 


from [1, N1]. After testing, Np was set to be approximately 30—100 depending on the size of the cluster. 
Fscale and Pcr were set to be 0.5 and 0.1, respectively. The convergence of the defects configuration with 
respect to r was checked by assigning r to be several values, e.g. 5, 10, 15 and 20 A. If the obtained 
configuration does not change with increasing r, a small r is accepted. To reduce the number of candidate 
states for a Vn, the shape of Regl was also considered according to the possible shape of the Vy. For 
example, for a bulk Va, Regl could be assigned to be spherical, while it is a rectangle at the GB given its 
preferred growth direction along the GB. 

(b). Mutation 

For each individual, a mutant vector (h’ PEE np ) Was generated according to (Fig. 1(c): 

h’ = round {x" + Fscale x(x" —x")}, (3) 
where random indexes r1, r2, r3 were chosen from [1, Np] and round is a rounding off function. Since the Va 
energy was optimized over integer space (Eq. 2),/ in Eq. 3 was rounded off. 

(c). Crossover 
To increase the diversity of the population, crossover was introduced [39]. A trial vector 
YS v, ..., v 


e. Va 


]was generated according to the following equation (Fig. 1(c)): 


(4) 


vV 


fe. hj if rand(0,1) < Per 
` |x}, if rand(0,1) > Per 


where i and j refers to the certain single defect in the cluster and individual in the population, respectively. 
Their values lie in the range ofie [1,n]and je [1, Np], and rand(0,1) is a uniform random number [0, 1]. 
(d). Selection 

Whether or not v’ replaces x/ as a new member of the population depends on the relative value of energy 
E(v’) and E(x’) (Fig. 1(c)). IfE(v’) is smaller than E(x’), x/is set tov’. Here, the energy of each offspring 
was defined to the system energy, which could be readily transformed into the corresponding Vn formation 
energy. After decoding a certain offspring, we obtained the corresponding unrelaxed V, configuration and 


created the V, in an atomic system with an approximate size.. 


To improve the efficiency for searching the ground state of the defects cluster, several methods were 
adopted. First, the individual with two identical components was avoided by means of penalty function. 
Equivalently, as the two Vs in a Vn occupy the same state, the corresponding individual was given a very 
high energy. Second, the relaxation of a V, expressed by an individual was performed every a specific steps. 
In principle, the system with a given initial rigid V, has to be relaxed at 0 K or at a finite temperature to get 
the accurate energy for the V,. Nevertheless, we found it is very time consuming as the system was relaxed 
at each step. To achieve a balance between the algorithm efficiency and accuracy, the system was only 
relaxed at a certain intervals (Fig. 1(d)), e.g. 20 steps in the present calculations. As relaxing the system, the 
conjugate gradient method was used to minimize the system energy at 0 K; the optimal time-step that leads 
to the maximal reduction in the system energy during each relaxation step was determined by using the 
golden section search method. In certain case, after relaxing the system with a Vm, a V in the Va 
spontaneously jumps to other neighboring sites. Correspondingly, the mapping relation of the relaxed space 
to the un-relaxed space was modified. As finding the stable V, configuration at the Fe GB, we found that the 
system has to be relaxed at a finite temperature, e.g. 300 K to obtain a low-energy V configuration. Third, to 
reduce the calculation time, the energies of all the individuals in the population were calculated in a parallel 
way. Finally, once the stable configuration Vno has been obtained, a large V, was optimized based on a small 
“core parts” Vno (Fig. 1(b)). In this case, only the configuration Vn-no was optimized by using the above DE 
algorithm (given the possible reconstruction of the Vn, nO=n-k with k>10 or other larger value). We also 
checked the parameter range of each individual and made periodical adjustment of the parameter to ensure 
that each component of the individual lies in the range of [1, N1]. We also tried to express the energy of a Vn 
by its geometrical parameters, e.g. the bonding between V-V and the bonding between the V and the 
surrounding atoms. Once the express was obtained, the time-consuming relaxation of the system was 


avoided. 


2.3. Energetics of Vn 


To characterize the energetics of a Vn, we defined the V, formation energy in a normal way as: 

Ey = Ey’ — Ey’ +nE gor» (5) 
where £y” is the energy of the system with a Vn, Z,” is the energy of a pure system without defects, E, is the 
cohesive energy per atom in the bulk of -4.122 eV. The clusters formation energy acts as an important 
parameter for determining the critical nucleation size of a cluster [26]. Particularly, AG? =—nk,T In S, + E; ; 
dAG? / dn =—k ,T In S, + AE," [26]. For a given temperature (T) and supersaturation concentration of the V 
(Sy), we can immediately obtain the critical nucleation size for a bulk void according to the energy level of a 
V near the void. 

To characterize the binding tendency of a V to the nano-void, we also defined the V-V, binding energy 
as: 
E; = Ep th, -Ej =E} AE, , (6) 
Considering that E/ is the energy level of the isolated V, the energy level of V near a cluster ( E; ) was 
thus defined as: 


Er =i, =F; (7) 


n+l 


2.4. The kinetics of a nano-void near the GB determined by using lattice kinetic Monte Carlo (LKMC) 
method 

Given the complex local atomic environment near the GB, the kinetics of a V in the vicinity of the GB 
was investigated by using LKMC method [42]. Contrary to the OKMC with a pre-constructed rate catalog 
for the events, in the present study, the rate catalog was rebuilt at every step of the LKMC. Each transition of 
the V, from the current state was assumed to involve only one jump of the V to its first- and/or 
second-nearest neighbor position. For each transition identified, the transition path was optimized by using 
the standard nudge elastic band (NEB) method [43], as illustrated in Fig. 1(e). The energy barriers of all the 
transition paths were calculated in a parallel way. For a large nano-void, this could greatly improve the 


calculation speed. For the transition with the energy barrier of Ea, its transition rate was calculated via 


E, ), where v, often in the range of 10’? —10" / s was assumed to be a fixed value of 10'°/s. The 


r =v, exp( 
B 


Boltzmann constant kg has a value of 8.617x10° eV/K. T denotes absolute temperature. After the rate 
catalog was constructed, the transition pathway was selected according to the standard KMC algorithm, 
where the chosen probability of a transition is proportional to the transition rate. Then, one V of Va jumps to 
the corresponding site. 

The diffusion coefficient (D) of a vacancy cluster at an elevated temperature was determined as the 
proportionally factor between the mean squared displacement of the mass center of the void and the time. 
We found that a large Vn is easily trapped in a local basin at a low temperature. Therefore, the D of a V, at a 
low temperature was extrapolated from the value at a high temperature once In(D) and the reciprocal of 


temperature follows the Arrhenius relationship (Fig. 1(f)). 


3. Results 
3.1. Formation, trap, leakage and emission of the V near Fe GBs 

We begin with investigating the behavior of single V near Fe GBs although the energetics of a V near Fe 
GBs has been studied systematically [47,48] and its kinetics near several GBs have also been studied 
[30,49-52]. By calculating the vacancy formation energy (Ey) at different sites near the GB, the site with 
the minimum Ey; (mE," ) was found, which acts as the stable position of a V at the GB. The V distributes 
within several layers near the GB (supplementary Fig. S1). Figure 4(a) presents mE, as a function of 
misorientation angle (@). Although there is no obvious trend inmEyy with respect toO mE; for some 
special low-X CSL-GBs (e.g. 25(2 1 0)/[0 0 1] and £5(3 1 0)/[0 0 1]) is noticeably larger than that for the 
GBs nearby. The twin boundary of £3(1 1 2)/[1 1 0] has the highest mE; . The low-angle GB with 0 less 
than 15° owes mE larger than that for the high-angle GB nearby, while the low-angle GB with @ higher 
than 75° has mE smaller than that for the high-angle GB nearby. To reveal the distribution feature of the 


stable site for the V at the GB, the spacing between these V-trapping sites (as marked by symbol d in 


supplementary Fig. S2(a)) is also presented in Fig. 2(a) as a function of @. It can be seen that d basically 
exhibits a U-shape curve with 0 , although it is scattered around a small value at the range of about 40—70° . 

In addition to the above static interaction parameters of the V with/within the GB, the LKMC and 
OKMC methods were further employed to explore the dynamic behavior of the V near the GB at an elevated 
temperature (mainly high temperature). Given the anisotropic strain field near the GB, lattice kinetic Monte 
Carlo (LKMC) was firstly employed to study the V trapping by the GB and reveal its dependence on the 
local atomic environment. It was found that the V interacted with different GBs via three distinct 
mechanisms. The V near X5(3 1 0)/[0 0 1] firstly followed a random walk and then moved straightforward 
towards the GB; it finally located at the region of the several atomic layers near the GB (Fig. 2(b)). Such 
behavior was readily understood based on the bulk-like migration energy curve in the region far away from 
the GB and a downhill shape of the curve near the GB (Fig. 2(b)). Near the low-angle GB of £85(7 6 0)/[0 0 
1] (0=81.2° , effective 0 =8.8° ), the V in a wider range underwent a non-random walk drift towards the GB. 
The longer-ranged strain field near the low-angle GB could give rise to such behavior, consistent with the 
finding in Ref. [53] that the low-angle copper GB has large sink strength due to the effect of GB stress. 
However, as the misorientation angle further decreases (e.g. £1741(30 29 0)/[0 0 1] with effective value of 
1.9°), the V near the region lying between the dislocation cores passed through the GB (Fig. 2(c)); similar 
phenomenon was observed for the V near £221(11 10 0)/[0 0 1] (supplementary Fig. S3). Such behavior was 
termed as the leakage effect. The V near the dislocation core of £1741(30 29 0)/[0 0 1] exhibited similar 
migration-trapping behavior to that near X85(7 6 0)/[0 0 1]. 

To measure the effect of the leakage on the V trapping, we further calculated the probability of a V 
moving away from the GB as a function of its distance normal to the GB at several GBs at 300 K (Fig. 2(d) 
and (e)). Initially, a V was put at the site in the bulk region about 25 A from the GB. As there is no a sink 
near the V, the V has the same probability moving forward and backward. As expected, near a sink, the 


probability for the V moving away from the GB is less than 0.5 due to the trapping of the V by the GB. The 


probability curves for the three GBs of ©5(2 1 0)/[0 0 1], X53 1 0)/[0 0 1] and £3(1 1 2)/[1 1 0] are similar. 
However, the probability in a wider range is less than 0.5 for the low-angle GBs of £85(7 6 0)/[0 0 1] and 
X221(11 10 0)/[0 0 1]. Yet, at an identical position near the GB, the moving-out probability for the V near 
the low-angle GBs (particularly near the GB of £1741(30 29 0))/[0 0 1]) is larger than that near the 
high-angle GBs of £5(2 1 0)/[0 0 1], £5(3 1 0)/[0 0 1] and £3(1 1 2)/[1 1 0]. 

To investigate the segregation of the V/SIA to the GB and their re-emission out of the GB, object kinetic 
Monte Carlo (OKMC) method was also employed. For this purpose, the atomic processes embedded in the 
OKMC model only include the V/SIA migration near the GB, trapping by the GB and emission from the GB. 
Figure 2(f) presents the V motion trajectories near several GBs at several temperatures of 300, 400 and 600 
K. At 300 K, the V was observed to be occasionally reemitted from the twin GB of £3(1 1 2)/[1 1 0] at 300 
K after segregation to the GB (Fig. 2(f)). As temperature increases to 400 K, the V was frequently remitted 
from the twin GB (Fig. 2(f)). At a high temperature of 600 K, the V trapped at the GB of &5(2 1 0) was 
found to be reemitted frequently out of the GB, while the emission is just an occasional event at the GB of 
X85(13 1 0) (Fig. 2(f)). The observed behavior is consistent with the magnitude of the V-GB binding energy 
in different GBs (the energy difference between Ey and Ey in Fig. 2(a)). Therefore, the single V behavior 
near the GB is related to both the local GB structure and the V-GB interaction parameters (e.g. V-GB 
binding energy). 

To further investigate the defect accumulation in NC/PC Fe, the stochastic cluster dynamics (SCD 
[54]//OKMC was used. The atomic processes within the GB were thus not allowed to occur, while in the 
grain interior all the atomic processes including diffusion, annihilation, binding, dissociation were allowed 
to occur. We calculated the trapping amount of the defects flowing into the GB. The typical values of V-GB 
binding energy (Eb or Ebjap) were adopted to represent several types of GBs. For the twin boundary of £3(1 
1 2)/[1 1 0], Ebiap for the V is 0.18 eV (Fig. 2(a)). For the special CSL-GB, Ebjap for the V in X5(2 1 0) of 


0.49 eV was chosen as a representative value. For general high-angle GB or low-angle GB, Ebja, for the V 


was given of 1.0 eV (the V-GB binding energy for £85(7 6 0)). As the leakage process was taken into 
consideration, the leakage rate was given of 0.5. Figure 2(g) presents the net trapping flux of the V/SIAs as a 
function of radiation dose under the condition of T=600 K, L=100 nm and r=107 dpaes”' . It can be seen that 
the accumulated V at the GB depend on Ebyap for the V. Compared to the other GBs, the net V concentration 
at the GB (C,”) of £3(1 1 2)/[1 1 0] is exceptionally lower. It implies that the twin boundary could not trap 
defects efficiently. For a low dose (e.g. 0.005 dpa), there is no significant difference in the defect 
accumulation in the other GBs due to the rarity of defects emission at a low dose. The C9” in the general 
GB/low angle GB is highest and becomes smaller as the V was leaked. TheCy”’ in the CSL-GB has an 
intermediate level in these systems. Therefore, the GB trapping efficiency for the V is sensitive to the V-GB 
binding energy. Consequently, as abundant Vs are present in the materials (e.g. under radiation), the quantity 


of Vs that segregate to the GB and compose the nano-void largely depends on the GB character. 


3.2. V-clustering at Fe GB at high temperature 

In addition to the relation of the trapping amount of the defects flowing into the GB to the local GB 
structure and V/SIA-GB binding energy, the determination of the defect accumulation at the GB also relies 
on the information of the defects within the GB, including the V kinetics at the GB, and V-V binding at the 
GB; of course, the V formation at the GB (Fig. 2(a)) also contributes the V accumulation within the GB (Fig. 
2(g)), although it also determines the V segregation via providing the V -GB binding energy (Fig. 2(a)). 

For the V kinetics at the GB, we calculated the V migration energy barrier at the GB by using the 
standard NEB method [43]. Figure 3(a) shows the V migration energy barrier (E4) along two nonequivalent 
paths within the GB as a function of the misorientation angle 0 . It suggests that E, along one path is much 
lower than that along the other one, as indicated by the V migration potential energy surfaces along the two 
types of paths at the GBs of £5(3 1 0) and £85(7 6 0) (supplementary Fig. S4). It implies that the V prefers 
to migrate on a certain low-energy barrier direction. Such speculation is supported by the quasi-one 


dimensional trajectory of the V motion at the GB drawn from the OKMC simulations of the V/SIA diffusion 


at the GB plane at 300 and 600 K (supplementary Fig. S5). Only at a high temperature, e.g. 600 K, the 
V/SIA at the special CSL-GBs of 25(2 1 0)/[0 0 1] and 25(3 1 0)/[0 0 1] and also X11(3 3 2)/[1 1 0] diffuses 
two dimensionally. In addition, E, for the V along the GB exhibits U-shape dependence on@. Basically, it 
could conclude that the V migration energy barrier is lower for the high-angle GBs compared to that for 
low-angle GBs. Different kinetics of the V at GBs may affect their clustering within the GB. 

For V-V clustering/binding at the GB, we calculated V-V binding energy (E; " ) (Eq. 6), which acts as 
an indicator for their clustering tendency. As shown in Fig. 3(b), Ey” does not heavily depend on@. 
Nonetheless, £; ” in most investigated GBs is positive and higher than that for a Fe bulk V-V pair except 
that in several GBs with the tilt axis of [1 1 0], indicating the attraction of Vs at the most GBs. Figure 3(c) 


presents the relation of EY 


Tto Ey . A positive correlation between the two variables was observed, 
indicating that the GB with high Ey has more stable V clusters within itself. The formation energy of small 
V clusters and the corresponding binding energy of the V with the V-cluster were also further calculated in 
several typical GBs (Fig. 3(d) and 3(e)). It could be seen that EY approximately show a linear relation to the 
cluster size of n, implying that the energy level of a V near the small V-cluster (the slope of the curve of the 
formation energy with respect to cluster size as defined by Eq. 7) basically remains a constant, different 
from the quick change of the V energetics near a bulk V, [20]. Furthermore, E; " exhibits a weak 
dependence on the cluster size of n, which arise from the linear configuration of the small V, within the GB. 
Nonetheless, the capillary law describes a different relation of the cluster energy to its size [54]. The V that 
repels to each other will be annihilated or locate at the GB randomly. In face-centered cubic (fcc) systems 
[55], the Vs randomly generated at GBs induced the GB phase transformation. In the present work, the 
emphasis is on the alternative path for Vs at the GB as formatting V clusters and GB voids. These energetic 
results guide the modeling of the V accumulation at different GBs. 


Given the net flux of the V into the GB and the Vs properties at the GB (Fig. 2(g)), the special CSL-GB 


acts as a candidate that probably develops a void structure under cumulative irradiation. Considering the V 


attracting each other (e.g. the GB of 25(3 1 0)/[0 0 1]), the trapped Vs could nucleate at the GB and grow to 
a GB void, as observed in experiments [3]. To explore the GB nano-void formation mechanism, the structure, 
morphology and energetics of nano-V clusters at Fe GB of ©5(3 1 0)/[0 0 1] were further explored by using 
the DE method, although the structures of small Vas (n<5) in bulk Fe are well-known [20,26]. Figure 3(f) 
presents the energetics of the V, at the GB. The energetics of the Vn have unique features. At the GB, from a 
single V to the 1D and then 2D configurations, the V has several typical energy levels corresponding to the 
specific folding configuration. The energy level oscillates with the size as n is larger than certain value, 
which is similar to the energetic behavior of small W bulk V,/SIA, in Ref. [56] (n<10 for Vn and n<20 for 
SIA, ). It arises from the formation of the interface between the defects and the atoms nearby. To understand 
the V, configuration transition during its growth at the GB, a mathematical setting for examining this 
transition was proposed (illustrated in supplementary Fig. S6(h)). The mode was based on the bonding 
among defects and the atoms nearby with a core idea to increase the shared area of the defect clusters and 
further increase the bonding strength. The formation energy for a defect cluster V2,/SIA2, with a 1D and 
folding configuration is E,=nE//,,-(2n-l)EX+E,, and E, =nE/!,,-[2(n—-NE} +(2n-lE?]+E,, , 
respectively. Here, E, / E,,is the energy increase in the configuration energy due to the formation of the 
interface between the cluster and the surrounding atoms; £} and Ef denotes the V-V binding energy in the 
two directions parallel and nonparallel to tilt axis, respectively. We assume that each isolated V/SIA in the 


two configurations has the same formation energy of E/!,,, Efra (Elsu = Ehu) and E =E. It is 


x 


energetically favorable for the cluster to be folded as E, = Z, , which leads to the condition 2n = = +1. For 


y 


b 


a different value of Ej/',,,, Efry we have the condition 2n > pp : 
E; + (Evisu ~ Evisa) / 2 

As displayed in Fig. 3(g), a single V resides at the GB in the form of V-I-V complexes (I denotes an 
interstitial) along [3 1 0] with the SIA located nearly on the GB plane and acting as the center of the two Vs 


sitting at the two symmetric sites on the two GB sides (Fig. 3(g)). Such complex has energy lower than the 


point-alike V structure found based on MS method [30,47—52]. Clusters V2—-V7 are composed of such V-I-V 
complex growing along [0 0 1] (Fig. 3(g)). Then, Vg is folded from 1D configuration along [0 0 1] to a 2D 
platelet configuration composed of the same complex (Fig. 3(g) for Vg); two columns of (V-I-V)4 complex 
are parallel to each other. During the above each stages of the V, growth, the configuration folding decreases 
the V, formation energy effectively (Fig. 3(f)). The mode of the V, growth along [0 0 1] and then 
undergoing folding repeats as n increases (Fig. 3(g)). As a result, Vn grows two dimensionally within the GB. 
Yet, the direction of the complex shifts from [3 1 0] to [1 1 1]. In addition, it was also found that the lattice 
surrounding the V, of the V-I-V complex is heavily strained approximately in the direction normal to the GB. 
With n increasing, more atoms contract towards the Va. The maximal displacement is as large as 0.7 A. 
Compared to the small strain around a bulk spherical void [57] (the present EAM potential predicts the 
displacement around a bulk V, is less than 0.1 A after relaxation), the strain here observed near a GB Vn is 
significantly large, which may arise from the fact that the compress of a bulk 3D void needs more work 
compared to the compress of a 2D GB void. A bulk void is a neutral sink because of a lack of associated 
strain fields [8,26]. Therefore, the Fe-GB nano-void could act as a preferential sink. The strain filed pattern 
also changes with n. The contraction strain field around the V, in the form of V-I-V complex is probably an 
intrinsic feature in the 2D defect GB system. It may be related to the elastic modulus. Since the greater the 
elastic modulus, the greater the stress needed for a given strain, the V, at a W CSL-GB may not have such 
low-energy V-I-V complex structure (the preliminary MS calculation of the two types of structures suggests 
that the V-I-V structure is unstable at W GB of £5(3 10 0)/[0 0 1)). 

As n further increases, the size of the cluster in the direction of [0 0 1] is so large that the DE fails to 
optimize the V, configuration as viewing the V, as a collection of single Vs. We turn to optimize the GB 
structures with different columns of Vs along [0 0 1] by the DE method. By such approximation, it implies 
an infinite size of the Va in the [0 0 1] direction, which may be rational only for a sufficiently large GB void. 


In this way, the optimized problem is finding the low energy configuration as creating Nc columns of Vs 


from Mc columns of atoms along [0 0 1] in the optimized region (equivalent number of Vs is Nx x Nc with 
Nx as the number atoms of a atom column along [0 0 1]). The optimization was carried out within a 
rectangle region with a maximal size of about 60 A in the direction of [1 3 0] and half-width of 10 A in the 
direction of [3 1 0] (the pink rectangle in supplementary Fig. S6(a)). It was found that for different values of 
Nc, the Vs at the GB developed to two types of structures. For a small Nc < Nc’ (or equivalently the void 
height A in supplementary Fig. S6(c) is smaller than a specific value), the vacant region is unstable. During 
the structural relaxation of a GB system with a void, a group of atoms surrounding the region rearranged 
their positions, accompanied by the contraction of the atoms in the grain interior towards the GB. As a result, 
the metastable GB phase was formed (supplementary Fig. S6(a) and S6(b)), which has a different structural 
unit to the pure GB. Actually, the sign of the transition of a large Vn to the metastable GB phase also found 
during optimizing the V, structure at the GB. For a large Nc > Nc’, the vacant region no longer collapsed, 
but formed a stable GB void phase with naked low-energy surfaces decorated with steps (supplementary Fig. 
S6(c) and S6(d)). Meanwhile, the strain near the GB recovered from the contraction strain (supplementary 
Fig. S6(e) and S6(f)). The void structure at the GB is similar to that in a bulk iron. 

We also found that the V-V binding energy at the GB and V kinetics at the GB are insufficient to 
determine the GB void morphology. For instance, the migration energy barrier of the single V along the GB 
and particularly V-V binding energy in some GBs (25(2 1 0)/[0 0 1], 25(3 1 0)/[0 0 1], 213(5 1 0)/[0 0 1], 
213(3 2 0)/[0 0 1], 225(4 3 0), 285(7 6 0), X85(13 1 0)/[0 0 1], 211(3 3 2)/[1 1 0]) have similar level. As the 
input parameters to the OKMC model were chosen from these GBs, we obtained similar V, 
nucleation/growth behavior for these GBs (supplementary Fig. S7). However, the DE calculation of the Vn 
configuration within theses GBs reveals two types of Va morphology at the GB. Only at the GBs of 25(3 1 0) 
and 25(2 1 0) does the V, have 1D-2D configuration, which further extends to the bulk region. The V, at the 
remaining GBs possesses a 1D configuration along the tilt axis except that has V-V repulsion interaction. 


These results suggest the important role of the GB character n coarse-grained modeling defect accumulation 


in PCs/NCs. Besides, the information of the large V, structure plays a role in predicting damage structure 


evolution at the GB in addition to the energetic and kinetic properties of point defects and small clusters. 


3.3. Structure and energetics of the W/Mo GB V, 

We have simulated V accumulation in Fe GBs (Figs. 2 and 3), finding that the behavior of clusters at the 
GB is essential to predicting radiation response of the GB. To explore the generalization of the results 
(mainly the V, formation at the GB and its relation to the local GB structure), we further calculated the Vn 
configuration at Mo/W GBs by using the DE method. Figure 4(a—d) presents the V, energetics at several 
Mo/W GBs (the V, configuration in bulk W is shown in supplementary Fig. S8; some of Vas in bulk Mo 
have identical configurations). It is found that the GBs could be divided into two categories according to 
their Ey’; as separated by the orange dash line in Fig. 4(a) and (c). The E at the GB of 25(3 1 0) and 25(2 1 
0)/[0 0 1] is significantly higher than that at the other GBs. The single V formation energies at the two GBs 
are also higher than that at the other GBs. Furthermore, the GB V, formation energy follows a different 
power law with respect to cluster size compared to that for a bulk V,. It is well-known that the formation 
energy for a bulk 3D V, has a dependence of n?’ [26] and the calculated values for bulk V, in W/Mo indeed 


1/2 


obey such relation. However, Ee was found to be expressed by the linear combination ofn'” andn', which 


ensures the convergence of Ey; > Ee |n > n0 (1/2<2/3<1), implying that a GB V, finally extends to the 
grain interior and becomes a bulk-alike void as its size is sufficiently large (supplementary Fig. S6(c) and 
S6(d)). Besides, the coefficient forn’? is much larger than that forn' (e.g. for Vn at W GB of 25(3 1 0), the 
two coefficients for n'? andn' are 4.21 and 0.46, respectively), suggesting that a small V, at the GB could be 
considered as a planar loop. The differentiation formation energy was clearly divided into two stages (Fig. 
4(b) and (d)). Both AES and AE decline quickly from 1 to about 5 and then oscillate between several 
typical levels. Such behavior could not be captured by the above n-power law continuum model which 


predicts that ABM = AE} [n —(n-1)""]/ [2 - 1] (note that in Ref. [20], 


E,(n) = E, +[E,(2)-£,|[n”” -(n-1)""]/[2*° -1]; here we have modified some symbols according to Eqs. 


(6) and (7)). In the following section, we will see that the V, with a low V energy level is more 
thermodynamically stable which serves as void nucleation center (the n that satisfies k,T7InS, = AE;"* 
determines the critical nucleation size [26]). 

Correspondingly, the V, configurations at the investigated nine GBs are clearly classified to two 
categories. In the W/Mo GBs of £5(3 1 0)/[0 0 1] and 25(2 1 0)/[0 0 ] (Fig. 4(e-h)), the V, initially adopts a 
linear configuration along the tilt axis on one side of the GB (note that the small V, at the two Fe GBs has 
a V-I-V complex configuration). The V, occupies the stable site for a single V at the GB, except for V2 at W 
GB of 25(2 1 0)/[0 0 1] (Fig. 4(f)). Then, the V, develops to a 3D structure centered at the GB plane, 
occupying additional sites in addition to the V stable site. During the development of a large Vn, it did not 
always grow on the basis of Vy.1, the reconstruction of the Vn was accompanied (e.g. from V4 to V5 in Fig. 
A(e) and for some bulk Vas in iron as suggested in supplementary Fig. S9). For some special clusters with a 
certain n, its configuration is mirror-symmetric about the GB plane, which may arise from the intrinsic 
symmetry of the CSL-GB. With the highly symmetric cluster as a center, the newly arriving Vs were 
absorbed at one side of the center until another symmetric V, was formed. Through such growth process, the 
Va extends in the three directions. We also found that the V, at the GB of 25(3 1 0)/[0 0 1] more frequently 
develops to a symmetric structure than the V, at the GB of 25(2 1 0)/[0 0 ]. In the remaining GBs, the Van 
was aligned along the tilt axis one-dimensionally with a V-I-V complex configuration. 

We then examined the relation of the GB V, morphology (3D versus 1D structure and point structure 
versus V-I-V complex) to the local GB structure. We note that a bulk V has equivalent neighboring sites with 
identical vacancy formation energy. Besides, the GB V, with cluster formation energy close to bulk values 
possesses a 3D configuration. Therefore, the configuration of the V, could be related to the fluctuation 
degree of the potential energy surface near the Vn. Figure 5(a) presents the V formation energy distribution 


at several GBs of W. Indeed, in the GBs of 25(3 1 0)/[0 0 1] and 25(2 1 0)/[0 0 1], the stable sites for the V 


with minimum formation energy are very near to each other not only in the direction along the tilt axis but 


also in the other direction. In the other GBs, the stable sites are aligned along the tilt axis. Furthermore, there 
are several atom columns with high V formation energy between the stable sites in the other direction along 
the GB. To further measure the flat degree of the potential energy surface at the GB, we calculated the 
standard error of vacancy formation energy (E/ ) for the group of GB atoms with low E/ as a function of the 
GB type in W. As shown in Fig. 5(b), the standard error well indicates two categories of the GBs with 
3D/1D V, morphology. A small error implies that the V prefers to condense along a certain direction, 
forming a 1D structure, while a large error implies that the local environment varies gently and the V, easily 
extends towards the surrounding sites via a 3D manner. 

In addition to the 3D/1D configuration of the V, at the GB, we further analyzed the strain around a GB 
V, and its relation to the GB structure (Fig. 5(c—e)). We found that the configuration of the V at the GB as a 
point or V-I-V complex is related to both the GB excess volume and the bulk modulus. The system either 
with low bulk modulus (e.g. Fe) or with small excess volume (e.g. W/Mo) has low energy V configuration in 
the form of V-I-V complex, while the V behaves like a point at the GB lattice site with large excess volume 


and high bulk modulus. 


3.4. Discrete model for the energetics of a GB V, 

As mentioned in the above section, the continuum model fails to capture the oscillated feature in the 
differential formation energy. Consequently, the thermal stability of a nano-void could not be well described 
by the model; the further formation and growth mechanism of the nano-void predicted by the model could 
be different to the real path. Next, we built discrete energetic models for GB V, by using the locally 
structural information of a cluster rather than only using the cluster size as in the continuum model. As 
observing the configuration of the V, (Fig. 6(a)), we found that the Vs within the V, were naturally 
connected by a series of bonds. Considering the interaction among Vs of the V, and interaction between the 


Vn with the surrounding atoms, the Vn formation energy was represented as the linear combination of the 


V-V and V-atom bonding: Ee = SO yvi tnya y-a) (Fig. 6(b)), where the sum extends to all the 


neighboring V-V and V-atom pairs;n,_,,is the number of i-th nearest V-V neighboring pair, /,,,1s the 
bonding strength for the V-V pair, n,_,,is the number of i-th nearest V-atom neighboring pair, and J;,_,, is 
the bonding strength for the V-atom pair. Note that, for a specific neighboring distance, Jj, and Jy.4; are 


constants. Such a linear model leads to the representation of the first-order differentiation of the formation 


energy as AE we = > (An, yf y-yi + Any_4Jy_4;)- For the bulk Vn, the V-V pair was restricted to the third 


i 


nearest neighbor, while V-A pair was only set to be the first nearest neighbor. For the GB Vj, the all the V-V 
and V-A neighboring distances were firstly calculated. According to their distribution, the typical values 
were chosen. In the discrete model, the parameters Iy-ņand Jy.4; could be obtained by fitting the calculated 
void formation energy to its structural variables ny.y; and ny-.4;. By considering the contribution of different 
neighboring V-V and V-atom pairs, the model naturally distinguishes the anisotropic feature of a nano-void 
with different facets. Besides, the model gives a discrete expression for the surface energy per area. The iron 
<1 0 0> surface is taken as an example. For this surface, E, , =n(2/,_,. +4J,_,,), while the surface area is 
S=ns, with sọ as the area of an atom. We could obtain the surface energy per area as 
y=E, -/S= (20, yy +4J5,_41)/ 5). Such model givesya value of 1.71J/m’, well consistent with the 
calculated value of 1.78 J/m? by using molecular statics. 

As presented in Fig. 6(c), such discrete bonding formula could describe both formation energy and its 
first-order differentiation with a very high precision (e.g. the maximum residual error for Fe bulk V, 
formation energy is as low as 0.05 eV in the discrete model based on V-V and V-atom bonding, which is as 
high as 1.0 eV in the continuum model). Consequently, the oscillation in the V energy level (Fig. 6(c)) was 
well reproduced by the discrete model, while the continuum model [54] neglecting the topology of the 
cluster and its structural relation with the underlying lattice fails in describing such oscillation feature (Fig. 
6(c)). In addition, the discrete model also naturally predicts a limited number of energy levels of the V near a 


void since the An,,_,,and An, ,, present in its energetic model have limited number of integer values. Besides, 


V—-Vi 


we could understand the discrete energy level of V based on such discrete model. A specific energy level in 


Fig. 6(c) corresponds to a local atomic environment that a new V occupies near a V,. Figure 6(e—g) suggests 
that the discrete model also applies to describe the energetics of the GB V,. Magana et al. [58] also proposed 
a modified Ising model to describe the clustering of Vs in bulk silicon voids. Although their model 
succeeded in describing the oscillation of V-voids binding energy, the model precision is lower than our 
model (e.g. the predicted binding energy for the V with V-V; is identical in Fig. 1 in Ref. [58]). Our 
discrete model also has predictive property. With the discrete model, the energy of the population in DE 
optimization of the V, structure (section 2.2) was directly calculated by only counting the number of 
neighboring V-V and V-atom pairs. Compared to the relaxation of the system with a Vj, the optimization 
speed was significantly improved. For a large V, (e.g. n>100), such discrete model based calculation of the 
V, energy could efficiently leads to the ground state of a void. 

The alternative discrete model was built based on the geometrical surface units of a Vn. A bulk Fe Vy 
contains two types of geometrical surface units with area of Ssa and s», respectively (as illustrated in Fig. 6(a) 


bulk 


by Vis). The Vn formation energy was thus represented by £,“ 


= Vays; , where n; is the number of surface 
i 


unit 7,7,1s its surface energy per area and s; is its area. But such model has a precision of only 0.7 eV to 
describe the bulk V, formation energy. The precision is too low to describe the differentiation formation 
energy that is on the order of 0-3 eV (Fig. 4(b) and (d)). Actually, a similar expression was employed in the 
Wull construction based energetic model for a void polyhedron [59], where the ratio of the surface energy 
pre area of different surfaces of the void determines the equilibrium shape of the void. The failing of these 
models to describe the energy of a nano-void arises from the neglected contribution of the junction between 
different surfaces to the total energy of the void. Compared to a very limited fraction of the junction in a 
meso-/macro-scopic void, the junction in a nano-void could play an important role in its energy. In a more 
approximated model based on spherical assumption of the shape of a void, the energy of a void is expressed 
as Ep“ = ys [26], wherey is its surface energy per area and s is its area. Such model neglects the 


anisotropic surface tension for different surfaces of the void in addition to the aforementioned the role of the 


junction in the total energy of the void. Therefore, these surface-based models do not apply to describe the 


formation energy of a nano-void, particularly its differential formation energy. 


3.5. Stability of a GB V, 

We finally investigated the thermodynamic and dynamic stability of a GB Va. A Va at a high 
temperature will dissociate. Correspondingly, the time for the V, to emit a V at an elevated temperature was 
used to characterize its thermodynamic stability (supplementary Fig. $10(a)). An emitted V from the GB Vn 
could locate at the GB or in the bulk (the corresponding time is t;/t2). Due to the low V energy level at the 
GB (Fig. 2(a), Fig. 4(b) and 4(d)) and migration energy barrier along the GB (Fig. 3(a)), a smaller time (¢7) 
was employed to characterize the lifetime of a V, at the interested GBs of 25(3 1 0) and 25(2 1 0) where the 
Va develops three dimensionally (Fig. 3(g) and Fig. 4(e-h)). As for the dynamic stability of the GB Vj, the 
diffusion coefficient of the Vn (D) within the GB was selected as an indicator. A small D implies that the Va 
remains immobile at a specific temperature. The D was calculated either directly via calculating the mean 
squared displacements of the V, at the GB at an elevated temperature or via the extrapolation from the 
high-temperature Arrhenius diffusion coefficient (D) (Fig. 1(f) and supplementary Fig. S10(b), section 2.4). 

Figure 7(a—c) shows the lifetime for a GB/bulk Vn in Fe/Mo/W. We found that the lifetime of a Vn 
depends heavily on the cluster size and its location (GB or bulk) in addition to temperature. It can be seen 
that although the lifetime for a V, exhibits an increasing trend with the cluster size, the increasing trend is 
non-monotonous. Some Vas with the special size (e.g. Ve, Vg and Vis, Vso in bulk Fe; V24, V31 at Mo GB of 
25(3 1 0)/[0 0 1]; Vi3 and V2; at W GB of 25(2 1 0)/[0 0 1]) have longer lifetime compared to the adjacent 
ones, acting as stable centers for void nucleation/growth. With the temperature increasing, the lifetime for 
both GB and bulk V, reduces greatly, indicating the instability of the V, at a high temperature 
thermodynamically. Yet, Vso in bulk Fe is still stable at a high temperature of 700 K, as indicated by a 
lifetime of over one day. The thermodynamic stability may be related to the structural symmetry of a Vn (Fig. 


A(e-h), Fig. 6(a)). At a low temperature of 300 K, most of Vas are stable except V; and V2. At an 


intermediate temperature of 500 K, only large cluster (n>25) and some small Vas (Vo, Vs, Viz, Via, Vis) are 
stable thermodynamically. We also found that the lifetime for a GB V, at a certain temperature is much 
smaller than that for a bulk Vn, implying that the bulk V, is much stable than a GB Va. In other words, the 
GB V, has to nucleate either on a large Vn or under high supersaturation of Vs at the GB, while a relatively 
small V, could at as the nucleation sites for a bulk V,,. Of course, the gaseous impurities also stabilize a void 
[7-9]. 

Next, the dynamic stability of the V, was obtained based on its D. Figure 7(d—f) displays the D for small 
Va at the GB/bulk in Fe/Mo/W. We found that the D for a Vn exhibits a decreasing trend with n, which 
however does not decrease monotonically with n and shows a discrete feature. The lower the temperature is, 
the more obvious the discretization in the D is (e.g. Ve and Vg have exceptionally small D at 300 K). The 
small V, still has a considerable mobility, as suggested by its moderate D. Similar relation of the cluster 
mobility to its structural symmetry was also found for bulk V, in iron (supplementary Fig. S11). The 
theoretical formula of a bubble D derived based on the surface diffusion [60] (the D is proportional to the 
reciprocal of R* or nf’, where R is the radius of a bubble) does not apply to a Vn, particularly at a low 
temperature. The D for a GB V, is not necessarily higher than that for a bulk Vn, although GBs have been 
considered as high-diffusivity paths [26]. For the GBs of 25(3 1 0) and 25(2 1 0) in Mo and GB of 25(2 1 0) 
in W, the GB Va has a higher D compared to that for a bulk Va, while the diffusion is slower along the Fe/W 
GB of 25(3 1 0) than in bulk Fe/W. 

Considering both the thermodynamic and dynamic stability, the stability of a V, is summarized, as 
illustrated in Fig. 7(f). At regime /, the V, is unstable both thermodynamically and kinetically (e.g. Vi—Vs at 
Mo GB of 25(3 1 0)). The V, with n<n0 falls into such regime, which contributes to the growth of a void via 
its dissociation (n>1) and diffusion. Such mechanism is different to the classic Ostwald ripening mechanism, 
where the large void grows only by the dissociation of small voids. Besides, the n0 for bulk V, is larger than 


that for a GB Van, implying that a GB void quickly becomes ripened compared to a bulk void. At regime 2, 


the Va with n=n*>n/ acts as a stable center for the V to nucleate both thermodynamically and dynamically 
(e.g. Vis and Vso in Fe bulk and V,3 and V2; at W GB of 25(2 1 0)). At regime 3, the V, (n0<n<n/) is 
thermodynamically stable but could diffuse (e.g. bulk Ve and Vs) to bind with other stable clusters. At 
regime 4, the V, is immobile but could dissociate, which could be the large cluster composed of some 


unstable branches and a stable center. 


4. Discussion 


4.1. Effect of the discretization nature of the energetics and kinetics on void formation and ripening 

We have shown that both the energetics and kinetics of the GB/bulk void exhibit a discretization 
feature (Fig. 3(e) and 3(f), Fig. 4(b) and 4(d), Fig. 6(c) and 6(e-g), Fig. 7, supplementary Figs. S10 and S11). 
Next we examine its effect on the void formation and ripening by using different binding models for the 
description of vacancy agglomeration. In the first case, we compared the model predicted average void size 
during SCD annealing to that detected in positron annihilation spectroscopy experiments on neutron 
irradiated single-crystal bcc Fe [61]. Results suggest that our discrete model leads to the average cluster size 
of 9 after annealing the system initially containing a concatenation of 10°% m° Vs at 50°C for 50 minutes 
(Fig. 8(a)), consistent with the experimental value, while the continuum model underestimates the average 
size. The further annealing simulations at high temperatures of 150°C and 250°C suggest that small clusters 
also drive the ripening process of the void in addition to free vacancies (Fig. 8(b) and 8(c)). For the GB void 
formation, there are no direct experiments as benchmarks to simulations. Yet, results indicate the significant 
differences in cluster size distribution predicted by the two different models (Fig. 8(d) and 8(e)), while the 
difference is small at a low temperature (Fig. 8(d)). The discrete model correctly reproduces the formation of 
particularly stable cluster centers, as predicted by the V energy level near a V, (Fig. 3(f), Fig. 4(b) and 4(d)) 
and also the cluster lifetime (Fig. 7(a-c)). Meanwhile, the maximum cluster size is smaller predicted by the 


continuum model than by the discrete model due to the underestimated energy level of the V near the Vn 


with low thermal stability (Fig. 3(e), Fig. 6(e-g)). 


4.2. Vacancy accumulation and void formation mechanism at GBs under radiation 

At a high temperature in a PC system, the GB response mechanism to irradiation depends on the GB 
structure [2—4,26,27]. The first determinate factor of the GB radiation response is the quantity of the Vs 
segregated to the GB, which is related to both the V-GB interaction strength and the distribution of the 
trapping sites at the GB (Fig. 2). Such factor dominates the vacancy concentration near the GB and further 
formation of the void-denuded zone near the GB. As illustrated in Fig. 9, the V energy level at the GB 
(E1—E4) varies significantly with the GB type (Fig. 2(a), Fig. 4(b) and (d)). These energy levels define two 


decisive variables of the V-GB binding energy (Eb 


trap 


) and the energetic driving force of the Vs to be 


clustered at the GB ( Eb 


clustering 


). Depending on the magnitude of Eb 


trap ? 


the V either permanently resides at 
the GB after segregation or is reemitted out of the GB (Fig. 9). The reemission of the trapped V at the GB is 
closely associated with the GB structure due to the scattered but relative small Ebja, for the V in different 
GBs (Fig. 2). The experimentally observed difference in the radiation response near different types of GBs, 
e.g. void-denuded zone formation and non-denuded zone mainly relies on the difference in the V-GB 
binding strength in these GBs. The twin boundary could not effectively trap the V due to the high V energy 


level (E/) and exceptionally small Eb 


trap 


(Fig. 2(a), Fig. 9); no sufficient defects could survive at the GB (Fig. 
2(g)). The low-angle GB and high-angle general GB could effectively absorb Vs from the grain interior (Fig. 
2(a), 2(e), 2(f)) due to the low energy level of the V at the GB (£3 and F4) and the corresponding large 
Eb, ay (Fig. 2(a)). The larger strain field near the low-angle GB also promotes the trapping of the V (Fig. 2(b) 
and 2(c)). Depending on the distribution of the trapping site along the GB, the trapping efficiency of the 
extremely low-angle GB is however lower compared to the general high-angle GB due to the leakage of the 
defects from the bulk-like region between the dislocation-centers (Fig. 2(d) and 2(g)). The special CSL-GB 
has an intermediate V energy level (Fig. 2(a), Fig. 9); correspondingly, the net quantity of Vs trapped at the 


GB lies between the values for the twin boundary and low-angle/general high-angle GBs (Fig. 2(g)). 


The additional determinate factor of the GB radiation response is their topology at the GB after the 
segregation. The relevant defects properties within the GB include migration style along the GB, migration 
energy barrier, binding energy between defects, and nucleation and growth processes of the Vn at the GB. 
The parameters for the nucleation (e.g. V-V, binding energy) are scattered with misorientation angle (Fig. 
3(b)), while the structure and kinetics of small Va depend on the GB structure (Figs. 3—5). The Vs collected 
at the GBs of 25(3 1 0) and 25 (2 1 0) develops to a compact 3D structure from a linear structure (Figs. 3 and 
3) via the migration of the single V (Fig. 3(a)) or/and small V, (Fig. 7(d—f)). The Vs at the other GB either 
could not bind with each other due to the repulsion between Vs or develop to a 1D structure lying within the 
GB with low thermodynamic stability (Fig. 3(b)). Besides, the growth of a GB void largely relies on Ep. A 


bulk V segregates to the GB energetically driven by the energy difference between Ey“ and E, (Eb 


trap ) = 


Once the segregated V contributes to the growth of the stable GB void, it is energetically driven by the 


i 


energy difference between Ey and Ey "ed (Eb ). Near a void of a certain size, L "remains a small 


clustering 
scattered value (Fig. 3(f)) and is close to zero averagely as predicted by the capillary law [26,54]. Therefore, 


Eb, and Eb 


trap clustering 


is mutually restricted by the equation of Eb,„ + Eb 


trap clustering 


< p Meanwhile, compared 
to other GBs, the special CSL-GB has a relative small Eb, (Fig. 2(a) and also as illustrated in Fig. 9). As a 
result, the defects accumulated at the low-angle GB and general high-angle GB with high Eb, would have a 


trap 


small energetic driving force of Eb, to grow (Fig. 9). In contrast, the defects collected at the special 


lustering 


CSL-GB have a large Eb 


clustering 


to further grow (Fig. 9). The Vs could finally agglomerate to the 3D GB void 


that is visible under TEM [3]. 


4.3. Effect of the coupling of grain size with GB character on vacancy accumulation at the GB 

The results in Figs. 2—8 are obtained in a single bi-crystal model. The real material contains the GB 
network composed of different types of GBs, connected by the triple junctions as observed in experiments 
[2,27,62]. Under certain kinetic conditions, radiation defects trapped at the GB could migrate along the GB 


network (the small migration energy barrier in Fig. 3(a) provides the kinetic condition), flowing about from 


one GB to another. Equivalently, the potential wells with different depths are connected (as illustrated in Fig. 
9 by the light yellow curve). The V coming from the grain interior initially occupy the wells with a high V 
energy level tend to move towards the well with a low V energy level (as illustrated in Fig. 9 by the yellow 
arrow). Correspondingly, the V accumulation in the material would be the collective consequence of 
multiple GBs. 

To support such speculation, we designed a system containing two sinks A and B with two energy levels 
of the V, as illustrated in supplementary Fig. S12(a). Sinks A and B could exchange defects with each other 
via emission, diffusion and trap processes in the OKMC model (supplementary Fig. S12(b)). The V was 
randomly produced in the two sinks at a rate. We found that the Vs in A with a high energy level flowed 
towards B before its reemission out of A. Moreover, the shorter the distance between A and B is, the more 
Vs finally flowed to B (supplementary Fig. $12(b)), implying the grain size effect in the process. Another 
SCD model where only the emission and trap processes were allowed to occur also confirmed the flow of Vs 
from A to B (supplementary Fig. S12(c)). Initially, the A and B were loaded with an equal number of Vs. As 
time passed, the Vs flowed from A to B with a lower energy level, as indicated by the decrease in the V 
percent in A and the increase in the V percent in B (supplementary Fig. S12(c)). Therefore, before the V 
could be emitted out of a local GB parts with small V-GB binding energy (or equivalently a high Ey ), the V 
finds its more stable site at the other GB (Fig. 2(a)). We termed this emission-inhibition behavior as the 
defects cruise within the GB, which significantly affects the V accumulation and distribution in a PC under 
irradiation. 

Next, the cruise condition was derived. As the time for V emission out of the GB 


(aha + 
e 


=f, ““s") ig longer than that for the defects moving along the GB to find a low energy 


( emission 


site in the other GB nearby (t,,,,. = L, / (2D) = L /(2D) = L / (2Dye *""“*") with L, as the distance that a 


cruise 


V has wandered during its diffusion process at the GB; here it was given of a typical value of grain size L), 


the V emission is inhibited. This leads to the following emission-inhibition condition: 


ERP +E (kyl Ep“ + Ep kyl . cae . ; : f i 
Per Em VED < 2D, tye‘ it Em MST) The physical insight manifestly encompassed in this equation is the 


coupling of the grain size (L) with the GB character (the characterization parameters EF and EF ), 


Vin 


Specifically, the cruise is easy to occur within the GB with a small Ef and E® for a given L or under the 


condition of a small L for a given GB with certain V properties (supplementary Table S1). For a large L and 


large Ei! and ES? 


Vm ? 


the emission tends to occur compared to the cruise (supplementary Table S1). For the other 


combinations of L, E and E” 


Vm ? 


it is not easy to arrive at a definitive conclusion on the relative rate of cruise 
and emission (supplementary Table S1). For instance, the V at £5(2 10)/[0 0 1] tends to cruise along the GB 
at 600 K as the grain size is less than 500 nm due to its small migration energy barrier of 0.36 eV, although 
its formation energy is as high as 1.22 eV (supplementary Table S1). 

The V cruise along the GB network would affect the relation of radiation-performance to the V-GB 
binding strength in PCs/NCs. The V emission plays a critical role in distinguishing the role of different GBs 
under irradiation by reducing the trapping efficiency of the GB for the defects in addition to the leakage (Fig. 
2(g)). As the V emission is inhibited under a certain condition, a definitive energy-space relation of the 
magnitude of the V-GB binding energy to the formation of the void-denuded zone near the corresponding 
GB (e.g. a large V-GB binding energy corresponds to the formation of the void-denuded zone) is broken. 
This brings about the uncertainty about the GB binding strength and radiation performance, as observed in 
NC Fe where different GBs of the same type exhibit difference in radiation performance at a high 
temperature [27]. Therefore, at a certain material parameter and radiation condition, the 
radiation-microstructure at a local region (particularly for a small grain size) is probably the result of 


collective behaviors of multiple GBs within a segment of the GB network. 


4.4. Implication of the V accumulation mechanism for designing the radiation resistance of a 
nano-structured metal 
Reducing the grain size is a well-known strategy for enhancing the radiation resistance of a poly-crystal 


metal. Our results provide new atomic insight into radiation damage reduction in NC metals from the view 


of the GB structure and have implications for designing and optimizing radiation-resistant NCs. The present 
work reveals several compromise approaches to tailor radiation performance of NCs and juggle the stability 


of the GB by taking advantage of the character of the GB itself. 


To improve the trapping efficiency of the defects by the GB, it qualitatively requires the increase in the 
binding energy of the defects with the GB. For example, the twin boundary binds with the V/SIA weakly 
compared to other GBs, and consequently the defect accumulation rate in the grain interior is larger than in 
other GBs at a high temperature due to the V/SIA emission from the GB. To obtain better radiation 
resistance of NCs at a high 7, the fraction of the twin boundary in the nano-structured system should be 
decreased. At a high radiation dose, the special CSL-GB is not recommended to have a large volume fraction, 
given the formation of the 3D void at the GB under irradiation (Figs. 3 and 4) because of the large energetic 
driving for the V to be clustered due to the relative high energy level of the V at the GB (Fig. 2(a)) and the 
flat vacancy potential energy surface at the GB (Fig. 5(a)). Besides, the void would further trap gaseous 
impurities, forming bubble and promoting boundary embattlement [8,63,64]. The general high-angle GB is 
recommended to give high volume fraction to improve the radiation-tolerance considering the difficulty in 


the GB voids formation. The large Eb 


trap 


gives rise to the high trapping capacity of the GB for the defects. 


Moreover, together with small Eb the trapped defects are hard to nucleate at the GB or difficult to 


clustering ? 
grow. Meanwhile, in view of the coupling of the grain size with the GB character (supplementary Fig. S12) 
and Table S1), the GB trapping efficiency is enhanced via the flow of the defects trapped at a CSL-GB with 
high defects formation energy to the general high-angle/low-angle GBs with low defects formation energy. 
Therefore, for high temperature application, the general GB with large defects-GB binding energy is 
suggested to have a large fraction together with a certain fraction of the CSL GB that stabilizes the GB 
network and meanwhile absorbs defects via defects transport towards the general GB, as pointed in Ref. [2]. 


There is a balance between the trapping and leakage of the defects near low-angle GBs. Although the 


low-angle GB shows good trapping ability for the defects (Fig. 2), the leakage effect emerges in the 


extremely low-angle GB, resulting in a low trapping efficiency in such system (Fig. 2(g)). Therefore, the 
system with a high fraction of the low-angle GB with appropriate angle should have high radiation 


resistance. 


5. Conclusions 

To sum up, based on the differential evolution, we developed a framework for determining the stable 
structure of a nano-void. Combining the framework with multiscale simulations, the structure, energetics 
and kinetics of both the point vacancy (V) and vacancy-cluster at iron (Fe), tungsten (W) and molybdenum 
(Mo) grain boundaries (GBs) have been systematically studied. We found that the void formation condition 
at GBs is closely related to the GB structure. At special coincidence-site-lattice (CSL) GBs of 25(3 1 0) and 
25(2 1 0), the V-GB binding energy is at a medium level, which although causes the V to be reemitted from 
the GB but also provides the driving force V-V binding at the GB. V-clusters with discrete size act as stable 
nucleation centers, developing from a linear configuration to a platelet and finally to a three-dimensional 
void with low thermodynamic stability. The growth is accompanied by the structural reconstruction of a 
group of vacancies. The small vacancy clusters with a relative mobility also contribute to the void formation 
in addition to the free Vs. The morphology of the V-clusters at the GB relates to both the excess volume and 
bulk modulus. General high-angle and low-angle GBs trap Vs efficiently (the V could be leaked from the 
exceptionally low-angle GB), where V-clusters only align one-dimensionally or hardly nucleate. A 
high-accuracy predictive linear model for the energetics of a nano-void both at the Fe/Mo/W GBs and in the 
grain interiors has been proposed based on the bonding among the vacancies and neighboring atoms of a 
nano-void. The model captures the anisotropic feature of a nano-void and reproduces the oscillation in the 
vacancy energy level near a nano-void, exhibiting superb advantages over widely used continuum model and 
Wulff construction based model for the energy of a meso-/maco-scopic void. Considering the flow of Vs 


from the GB with high V formation energy towards the GB with low formation energy, the collective 


behavior of GB networks plays a role in GB void formation. To render the poly-crystal Fe/W/Mo more 
radiation-resistant at a high temperature, the general GB with large defects-GB binding energy is suggested 
to have a large fraction. The present work provides new knowledge in understanding the GB void formation 
and enables sound basis for tailoring poly-crystal’s radiation resistance to void formation by reducing the 


fraction of special CSL-GBs. 
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R trapping sites at the GB. (b) The typical V motion trajectories near the GBs of £5(3 1 0)/[0 0 1] and £85(7 6 
x< 0)/[0 0 1] at 300 K drawn from lattice kinetic Monte Carlo (LKMC) simulation of the V diffusion near the 
i - GB. The horizontal line indicates the GB position. Axis Z is normal to the GB plane, while axis Y is parallel 

. to the GB plane. For 25(3 1 0)/[0 0 1], Y and Z is along [1 3 0] and [3 1 0], respectively, while for &85(7 6 
0)/[0 0 1] Y and Z is along [67 0] and [7 6 0], respectively. (c) Several motion trajectories of the V near a 
low-angle GB of £1741(30 29 0)/[0 0 1] at 300K. Here, symbols I1, I2 and I3 mark three initial sites of the 
introduced V, while El and E2 mark the final trapping sites. (d) and (e) The probability of a V moving away 
from the GB as a function of its distance normal to the GB at several GBs at 300 K. Initially, a V was put at 
the site about 25 A from the GB. (f) The motion trajectories of a V near the GBs at 300, 400 and 600K 
drawn from object kinetic Monte Carlo (OKMC) simulations of the V behavior near the GB. Initially, the V 


is about 6 nm away from the GBs. The horizontal line indicates the GB position. (g) The V concentration at 
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the GB (C ) as a function of radiation dose at 600 K for typical values of the V-GB binding energy (Eb or 


Ebirap), as marked in (a). Here, ratejea, denotes the leakage rate of the V from the GB. 
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Fig. 3. (a) The V migration energy barriers (E4) along two paths within the Fe GBs. (b) The V-V binding 


energy (E, ) at the GB as a function of the misorientation angle. (c) The relation of the V-V binding energy 


to the minimum V formation energy at the GB. (d) The V-cluster formation energy (E; oe z) as a function of 


the cluster size n. (e) The binding energy of the V-V, (E, ’”) at the GBs as a function of n. (f) The Vn 


formation energy (E;% ) (right axis) and its one-order difference (£;” 


-E;',) (left axis) as a function of n 


within the GB. Here, nl” the critical cluster sizes for the V, structural transition from a 1D configuration 


aligned along [0 0 1] tilt axis to a 2D configuration with two columns of Vs parallel to the tilt axis. Similarly, 


n2 is the critical cluster size for the V, transition from a two-column configuration to a three-column 


configuration along the tilt axis. Here, the black line indicates the result of Es; , predicted by the continuum 


model. (g) The ground state structures of several Vs within &5(3 1 0)/[0 0 1]. Here, dz is atomic 


displacement along [3 1 0]. 
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Fig. 4. (a) and (b) The formation energy for the vacancy cluster (Ey ) and the differential energy (AE; = 
Ej 7 =E ) as a function of n for several W GBs. (c) and (d) The E and AE; as a function of n for several 
Mo GBs. For comparison, the energetic results for bulk V, in W/Mo are also shown. (e) and (f) The 
groundstate configurations of typical Vas in W GBs of 25(3 1 0)/[0 0 1] and 25(2 1 0)/[0 0 1], respectively. 
(g) and (h) The groundstate configurations of typical Vas in Mo GBs of 25(3 1 0)/[0 0 1] and 25(2 1 0)/[0 0 


1], respectively. In (e-h), the V in a V, is represented by a cyan sphere. The first and second nearest V of a V 


is connected by a green and pink line, respectively. 
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Fig. 5. (a) The distribution of vacancy formation energy within several W GBs. (b) Standard error of 
vacancy formation energy ( E/ ) for the group of GB atoms with low E/ as a function of the GB type in W. 
(c) and (d) The structure of a three-dimensional (3D) and 1D vacancy cluster at the W GB of 25(3 1 0)/[0 0 
1] and 213(3 2 0)/[0 0 1], respectively. (e) Relation of the V, configuration within the investigated Fe/Mo/W 


GBs with the GB excess volume and the bulk modulus of Fe/Mo/W. 
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Fig. 6. (a) The groundstate configurations of typical Vas in Fe bulk. Here, the V in a V, is represented by a 
cyan sphere. The first and second nearest V of a V is connected by a green and pink line, respectively, while 
; dl and d2 marked near V4 denotes the respective neighboring distance. (b) The local atomic environment for 
a Vn described by its neighboring atom (Ai), neighboring vacancy (Vi). Here, n,_4; nyy denotes the 
number of V-Ai, V-Vi, respectively. q(n,_,,,%,_,,)andq(An,_,,,An,_,,) denotes the linear combination of 
Ny_4io Ny_y, and An,_,,,An,_,,, respectively. The green cubic and blue sphere represents a vacancy and its 


neighboring atom. (c) The formation energy for the vacancy cluster (EE ) and the differential energy (AES = 


Ef —Ef ) as a function of n for bulk Fe. (d) The relation of the V energy level to the local atomic 


n+l 


environment. The first, second and third nearest neighbor of the V is respectively connected by the red, 
green and blue line, respectively. The green cubic indicates the V added on a Vn. (e-g) Model prediction of 


E and AE; in W bulk, and GBs of 25(2 1 0)/[0 0 1] and 25(3 1 0)/[0 0 1], respectively. 
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Fig. 7. (a-c)The lifetime of a V, at several typical temperatures as a function of n. (d-f) Diffusion 
coefficients for small Va within Fe/Mo/W GBs of &5(3 1 0)/[0 0 1] and £5(2 1 0)/[0 0 1]. For compassion, 
bulk values for V, are also presented. In the inset figure of (f), E and K denotes the thermal and kinetic 


stability, respectively of a Vp. 
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=- Fig. 8. (a) The average size of the V» during annealing a Fe bulk system that initially contains a 
| concentration of 5x 10%m™ of single vacancies for 50 minutes at temperature of 50°C. (b) and (c) The 
distribution of the size of V, entering into the clustering and dissociation reaction at 150 and 250°C, 
respectively. (d) and (e) The cluster density distribution after annealing a system that initially contains a 
= concentration of 5 x10'*m”~ of single vacancies for 50 minutes for Fe GB of £5(3 1 0)/[0 0 1] and W GB of 


x5(2 1 0)/[0 0 1], respectively. 
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_ Fig. 9. Illustration of the defects trapping and clustering at a single GB and a GB network. The symbols 


E0-E5 denote the V energy level at different GBs. bya, is for the V-GB binding energy. Eb/ is for the V-V 


) binding energy at the GB, while £b2 is for the binding energy of the V with a GB void. The inequality of 


Eb. + Eb 


trap clustering 


< E;;" describes the restricted relation of the V-GB binding energy to the energetic driving 


force for the V-clustering at the GB. The gray bow-shaped curve around a void indicates the compressive 


strain near the GB void. The inequality of Pe” 


Wf +Eyn )/(kpT) 


(Ey! +Ern Mp 


< 2D e 


” describes the coupling 


condition of the grain size (L) and in-boundary V formation and migration when the emission of the V is 


inhibited due to the cruise of the V along the GB to find a low energy site. Here, £$? and E8? denotes the 


of Vin 


vacancy formation energy and migration energy barrier at the GB, respectively, while Ey" and Epp“ is the 


corresponding value for a bulk V. D, and tis the pre-exponential factor for vacancy diffusion coefficient and 


emission time at the GB, respectively. Here, D,t, =d,, withd, as jump distance of the V at the GB. T is 


absolute temperature and kg is Boltzmann constant. 


